#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

#include "Matrix.h"

#ifndef DDMCMC_Vector
#define DDMCMC_Vector

template<class Tp>
class Vector : public Matrix<Tp>{
public:
	// construction & destruction
	Vector(){};
	Vector(const Vector<Tp>& a);
	Vector(const int row, const int col);
	Vector(const int size){ SetRowVector(size); };
	~Vector(){};

	// operators
	Vector<Tp>& operator=(const Vector<Tp>& x);
	Tp&				operator()(const long index);
	Tp&				operator()(const int row, const int col);

	Vector operator+(const Vector<Tp> &x);
	Vector operator-(const Vector<Tp> &x);
	Vector operator*(const Matrix<Tp> &x);
	Vector operator*(const Vector<Tp> &x);

	friend std::ostream &operator<< (std::ostream &strm, Vector& vec){ return strm << vec.data << endl; };
	void operator>>(std::ofstream &strm);
	void operator<<(std::ifstream &strm);

	Vector		operator+(const Tp y);
	Vector		operator-(const Tp y);
	Vector		operator*(const Tp y);
	Vector		operator/(const Tp y);

	// operations
	void	SetRowVector(const int size);
	void	SetColVector(const int size);
	bool	RowVector(void);
	bool	ColVector(void);
	bool	ExpandTo(const int size);


	double	Norm(const double l);
	Tp		Projection(const Vector &a);
	Tp		DotProduct(const Vector &a);

	void	CopyContent(const Vector<Tp> &b, const int num);
	virtual Vector  Conv1D(const Vector& kern, const bool bcircular);	// circular convolution in 1D
};


////////////////////////////////////////////
// construction & destruction
template<class Tp>
Vector<Tp>::Vector(const Vector<Tp>& a){
	(*this) = a;
}

template<class Tp>
Vector<Tp>::Vector(const int row, const int col){
	SetDimension(row, col);
}

////////////////////////////////////////////
// operators

template<class Tp>
Vector<Tp>& Vector<Tp>::operator=(const Vector<Tp>& x){
	//beware of self assignments : x = x
	if (this != &x){
		if (rows() != x.rows() || cols() != x.cols()){
			SetDimension(x.rows(), x.cols());
		}

		data = x.data.clone();
	}
	return *this;
}

template<class Tp>
void Vector<Tp>::operator>>(ofstream &strm)
{
	// dimensions
	strm.write((char*)&nrow, sizeof(nrow));
	strm.write((char*)&ncol, sizeof(ncol));

	// data
	strm.write((char*)data.data, sizeof(Tp)*size());
}

template<class Tp>
void Vector<Tp>::operator<<(ifstream &strm)
{
	int row_read, col_read;

	MatrixFree();

	// dimensions
	strm.read((char*)&row_read, sizeof(nrow));
	strm.read((char*)&col_read, sizeof(ncol));

	if (row_read>0 && col_read>0)
		SetDimension(row_read, col_read);
	// data
	strm.read((char*)data.data, sizeof(Tp)*size());
}

template<class Tp>
Vector<Tp> Vector<Tp>::operator+ (const Vector<Tp>& x){
	Vector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary matrix
	temp.data = data + x.data;

	return temp; //return the new matrix
}

template<class Tp>
Vector<Tp> Vector<Tp>::operator+ (const Tp y){
	Vector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary matrix
	temp.data = data + y;

	return temp; //return the matrix
}

template<class Tp>
Vector<Tp> Vector<Tp>::operator - (const Vector<Tp>& x){
	Vector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary Vector
	temp.data = data - x.data;

	return temp; //return the new Vector
}

template<class Tp>
Vector<Tp> Vector<Tp>::operator - (const Tp y){
	Vector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary Vector
	temp.data = data - y;

	return temp; //return the Vector
}

template<class Tp>
Vector<Tp> Vector<Tp>::operator * (const Vector<Tp>& x){
	Vector<Tp>	temp(nrow, x.ncol); //intialize the temporary Vector

	cv::Mat a = data.clone(), b = x.data.clone();
	a.convertTo(a, CV_64FC1);
	b.convertTo(b, CV_64FC1);
	cv::Mat c = a*b;
	c.convertTo(c, cv::DataType<Tp>::type);
	temp.data = c;

	return temp; //return the result of multiplication
}

template<class Tp>
Vector<Tp> Vector<Tp>::operator * (const Matrix<Tp>& x){
	Vector<Tp>	temp;

	temp.SetDimension(rows(), x.cols()); //intialize the temporary Vector
	cv::Mat a = data.clone(), b = x.data.clone();
	a.convertTo(a, CV_64FC1);
	b.convertTo(b, CV_64FC1);
	cv::Mat c = a*b;
	c.convertTo(c, cv::DataType<Tp>::type);
	temp.data = c;

	return temp; //return the result of multiplication
}

template<class Tp>
Vector<Tp> Vector<Tp>::operator * (const Tp y){
	Vector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary Vector
	temp.data = data*y;

	return temp; //return the Vector
}


template<class Tp>
Vector<Tp> Vector<Tp>::operator / (const Tp y)
{  //multiply by a scalar constant
	Vector<Tp> temp;

	temp.SetDimension(rows(), cols()); //intialize the temporary Vector
	temp.data = data/y;

	return temp; //return the Vector
}

template<class Tp>
Tp& Vector<Tp>::operator ()(const int row, const int col){
	return data.at<Tp>(row, col);
}

template<class Tp>
Tp& Vector<Tp>::operator ()(const long index){
	return data.at<Tp>(index / ncol, index % ncol);
}

template<class Tp>
void Vector<Tp>::SetRowVector(const int size)
{
	SetDimension(1, size);
}

template<class Tp>
void Vector<Tp>::SetColVector(const int size)
{
	SetDimension(size, 1);
}

template<class Tp>
bool Vector<Tp>::ExpandTo(const int size)
{
	if (RowVector())
		return Matrix<Tp>::ExpandTo(1, size);
	else
		return Matrix<Tp>::ExpandTo(size, 1);
}

template<class Tp>
bool Vector<Tp>::RowVector(void)
{
	if (cols()>0 && rows() == 1)
		return true;
	else
		return false;
}

template<class Tp>
bool Vector<Tp>::ColVector(void)
{
	if (cols() == 1 && rows()>0)
		return true;
	else
		return false;
}

template<class Tp>
Tp Vector<Tp>::Projection(const Vector<Tp> &a){
	Tp sum = 0;
	for (int i = 0; i<size(); i++)
		sum = sum + Data(i)*a.Data(i);
	return sum;
}

template<class Tp>
Tp Vector<Tp>::DotProduct(const Vector<Tp> &a){
	Tp sum = 0;
	for (int i = 0; i<size(); i++)
		sum = sum + Data(i)*a.Data(i);
	return sum;
}

template<class Tp>
double Vector<Tp>::Norm(const double l)
{
	if (l == 1){
		return cv::norm(data, cv::NORM_L1);
	}else if (l == 2){
		return cv::norm(data, cv::NORM_L2);
	}else{
		int		i;
		double	sum = 0.0;

		for (i = 0; i<size(); i++)
			sum = sum + pow(fabs((double)Data(i)), l);
		sum = pow(sum, 1.0 / l);
		return sum;
	}
}

template<class Tp>
void Vector<Tp>::CopyContent(const Vector<Tp> &b, const int num)
{
	for (int i = 0; i<num; i++)
		(*this)(i) = b.Data(i);
}

// circular convolution in 1D
template<class Tp>
Vector<Tp> Vector<Tp>::Conv1D(const Vector& kern, const bool bcircular)
{
	Vector<Tp>	I;
	long		i, k, m, nsize, nsize_kernel;

	I.SetDimension(rows(), cols());
	I.Zero();

	nsize = size();
	nsize_kernel = kern.size();

	if (bcircular)
	{
		// circular convolution
		for (i = 0; i<nsize; i++)
		{
			for (k = 0; k<nsize_kernel; k++)
			{
				m = (i - (k - (nsize_kernel - 1) / 2) + nsize * 100) % nsize;

				I(i) = I(i) + Data(m)*kern.Data(k);
			}
		}
	}
	else
	{
		for (i = 0; i<nsize; i++)
		{
			for (k = 0; k<nsize_kernel; k++)
			{
				m = i - (k - (nsize_kernel - 1) / 2);

				if (m >= 0 && m<nsize)
					I(i) = I(i) + Data(m)*kern.Data(k);
			}
		}
	}

	return I;
}






#endif